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SUMMARY 


A  set  of  algorithms  has  been  developed  for  the  estimation  of 
strain-time  histories  in  1-D  spherical  fields.  Methods  for  the 
estimation  of  stress-difference  and  stress-bounds  have  also  been 
developed.  Qualitative  statements  regarding  the  effects  of 
possible  errors  in  method  and  instrument  measurement  error  have 
been  made  and  substantiated. 


Further : 


At  the  current  state-of-the-art,  Lagrangian  analysis 
must  be  considered  as  an  adjunct  to  laboratory  methods 
and  parametric  calculations  in  the  development  of  better 
material  models. 

The  methods  outlined  in  this  paper  require  no  explicit 
shape  function  to  be  assumed.  Thus,  complex  flow  fields 
can  be  dealt  with. 

Accuracy  of  strain-time  histories  are  affected  to  first 
order  by  errors  in  late-time  displacement  (as  due,  in 
practice,  to  baseline  shifts,  instrument  calibration, 
and/or  accuracy) . 

Although  only  briefly  discussed  here,  we  reiterate  a 
known  result:  stress-difference  time  history  estimation 
is  very  sensitive  to  errors  in  acceleration  and  dar/dh. 

Local  estimates  of  time  of  arrival  and  flow  s':-‘face 
should  be  used  in  estimating  strains,  unless  errors  in 
gage  surveying  or  gage  availability  dictate  otherwise. 
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A  more  global  estimate  of  time  of  arrival  and  flow 
surface  should  be  used  in  estimating  stress-bounds  in 
order  to  avoid  extrapolation. 

Practically,  difficulty  still  remains  in  validating  the 
degree  of  one-dimensionality  of  the  in  situ  environment. 
Because  the  state  of  current  instrumentation  does  not 
allow  a  measure  of  the  uncertainty,  the  only  clear  (but 
by  no  means  unequivocal  or  definitive)  method  would  seem 
to  be  analysis  of  a  corresponding  two-  or  three- 
dimensional  calculation. 

Difficulties  and  ambiguities  in  validating  measurements 
from  current  motion  and  stress  must  be  dealt  with  in 
order  to  remove  possible  first  order  errors  from  strain 
estimates . 


CONVERSION  FACTORS,  NON-SI  TO  SI  (METRIC) 
UNITS  OF  MEASUREMENT 


Non-SI  units  of  measurement  used  in  this  report  can  be  converted 
to  SI  (metric)  units  as  follows: 


Multiply 

By 

To  Obtain 

degrees  (angle) 

0 .01745329 

radians 

feet 

0 . 3048 

meters 

gallons  (US  liquid) 

3.785412 

cubic  decimeters 
( liters ) 

inches 

2.54 

centimeters 

kips  (force) 

4.448222 

kilonewtons 

kips  ( force  per 
square  inch) 

6.894757 

megapascals 

megatons  (nuclear 
equivalent  of  TNT) 

4  . 184 

peta joules 

pounds  (force)  per 
square  inch 

6.894757 

kilopascals 

pounds  (mass) 

0.4535924 

kilograms 

pounds  (mass)  per 
cubic  foot 

16.01846 

kilograms  per 
cubic  meter 

kbars 

98 

megapascals 

feet 


1/.3048 


meters 
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SECTION  1 
INTRODUCTION 


Estimation  of  stress  and  strain  paths  followed  in  geological 
media  using  measured  velocity  and  stress  gages  obtained  from 
purely  one-dimensional  events,  a  technique  commonly  called 
Lagrangian  Analysis  of  Stress  and  Strain  (LASS),  has  been 
employed  for  at  least  15  years  (Refs.  1,  2,  3,  4,  5,  6). 
Currently,  such  techniques  are  being  used  to  provide  insight  into 
the  information  needed  to  construct,  ideally,  more  precise 
material  models  for  use  in  2-D  finite  difference  calculations. 

The  results  from  recent  LASS  analysis  of  dry  soils  to  date 
have  not  been  self-consistent  nor  yielded  models  which  have  been 
shown  to  be  better  than  those  currently  arrived  at  by  laboratory 
tests  and  parametric  calculations.  The  reason  for  this  is 
twofold : 

First,  it  is  due  to  the  inability  to  validate  motion  and 
stress  measurements  to  the  required  degree  of  precision  required, 
resulting  from: 

1.  Lack  of  repeatability  for  similarly  emplaced  gages. 

2.  Imprecisely  defined  arrival  and  rise  times. 

3.  Inadequate  calibration  techniques. 

4.  Inability  to  objectively  correct  for  gage  drift  after 
shock  passage. 

5.  Inability  to  measure  "hoop"  stresses. 

6.  Limited  survivability  hampering  validation. 
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Limited  validation  of  geometrical  symmetries  required 
for  analysis. 
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8.  Inconsistencies  between  motion  and  stress  measurements, 
partly  due  to  inaccurate  early-time  response  of  stress 
gages . 

Second,  it  is  because  very  little  has  been  done  regarding 
the  influences  of  possible  systematic  and  inherent  errors  and 
errors  due  to  varying  judgment  as  to  what  constitutes  a  "good" 
data  set  and  appropriate  methods  of  data  reduction.  For  example, 
separate  LASS  analyses  of  gorund  motion  data  from  event  MP2 
(Refs.  7  &  8 ) ,  by  Trulio,  (Ref.  9),  and  Allen  (Ref.  10)  came  to 
remarkably  differing  conclusions  regarding  suitability  of  current 
models . 

This  paper  develops  another  method  of  Lagrangian  analysis, 
which  deviates  somewhat  from  earlier  "path-line"  (Refs.  1,  2,  3, 
4,  5,  6)  and  "traveling  wave"  formulations.  The  method  may  be 
considered  as  computationally  simplifying  the  path  line  method 
(Refs.  1  and  3)  in  that  directional  derivatives  along  an 
arbitrary  path  are  replaced  by  partial  derivatives  related  to  the 
time-of-arrival .  It  may  also  be  considered  a  generalization  of 
the  traveling  wave  method  (Refs.  4  and  5),  applicable  for 
arbitrarily  complicated  one-dimensional  flow. 

The  basis  for  the  method  is  a  transformation  which  aligns 
all  gage  records  at  zero  time,  an  idea  originally  suggested  by 
Seamen  (Ref.  3).  After  necessary  computation  of  the  required 
derivatives  at  a  point  on  the  appropriate  surface,  the  results 
are  transformed  back  into  original  gage-time  using  an  appropriate 
fit  to  the  time-of-arrival.  In  the  situations  where  evaluation 
of  a  definite  integral  along  a  functional  surface  is  needed,  the 
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integration  is  calculated  as  a  line  integral  on  the  transformed 
surface . 

Differences  in  interpolation,  smoothing,  curve  and  surface 
fitting,  and  solution  methodology  do  lead  to  markedly  different 
results  in  many  cases;  the  algorithms  that  were  finally  chosen 
are  explicitly  derived  and  justified. 

While  this  paper  and  the  subsequent  report  (Ref.  11),  which 
applies  the  techniques  herein  to  a  specific  event,  do  not  address 
all  the  possible  sources  of  inconsistency ,  it  will  point  out  many 
which  have  not  received  attention  before.  Specifically, 

1.  Errors  due  to  instrumentation  inaccuracy. 

2.  Methodological  errors  in  estimation  of  strains  and 
stress  bounds. 


The  equations  are  general  for  nonisentropic  waves2  and  are 
valid  under  the  following  assumptions  (Refs.  1  and  12): 

1.  Heat  flow  is  to  be  neglected. 

2.  All  body  forces  are  to  be  neglected. 

3.  No  energy  sinks  or  sources  are  postulated. 

4.  All  flow  is  spherical. 

Assumptions  1  and  2  are  considered  valid  because  of  the  very 
short  time  intervals  and  high  stresses  involved.  Assumption  3  is 
generally  valid  until  the  point  where  discontinuous  phenomena 
occur,  such  as  cracking  of  material  or  onset  of  phase  changes. 


Equation  1  expresses  conservation  of  mass.  Equation  2 
expresses  conservation  of  momentum.  Equation  3  is  the 
kinematic  particle  velocity.  Note:  the  convention  of 
Equation  2  is  that  a  tensile  stress  is  negative. 

A  rapidly  attenuating  shock  induced  by  a  detonation  on 
or  in  the  ground  being  the  case  of  present  interest. 
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Assumption  4  is  ideally  valid  for  a  perfectly  spherical 
detonation  in  an  isotropic,  homogeneous,  full  space  and  implies 
an  absence  of  tangential  motions  and  any  shear  stresses  due  to 
vorticity.  Perfectly  spherical  flow  can  be  interpreted  as  a 
singularity  at  the  origin  in  the  equations  of  motion;  i.e.,  as  a 
mass  source,  In  the  case  of  a  contained  detonation  at  large 
radii,  the  same  idealization  applies. 

Practically,  for  detonations  in  near  surface  layered  media 
with  measurements  of  unknown  absolute  accuracy  relatively  near  an 
expanding  crater,  the  posed  problem  is  really  a  complex  boundary 
value  problem  in  a  nonideal  medium.  Investigation  into 
disturbing  effects  resulting  from  these  two  considerations  and 
their  implications  for  solution  of  1-3  has  just  begun. 


SECTION  3 

APPROXIMATION  OF  FLOW  SURFACE 


Given  an  event  meeting  the  assumptions  in  Section  I  and  a 
radial  stress  and  motion  field  approximately  described  by 
observed  Lagrangian  gages,  Equations  1-3  can  be  inverted  for  the 
field  strains  and  field  stress  difference. 

In  order  to  do  this  we  consider  the  spherical  one¬ 
dimensional  flow  surface  in  Figure  la  approximately  represented 
by  the  observed  particle  velocity  variable3  (i.e.,  gage 
response),  u  as  a  function4  of  the  coordinates  h  and  t;  that  is, 

u  =  u ( h, t ) . 

Further 

u(ho,t:)  =  ft)  and  r(ho't)  "  ho  +  J  u(hG/t')dt'.  (4) 

h  t 

o 

Define  h  as  the  initial  gage  position,  and  t  =  f  (h) ,  as 

the  time-of-arrival  of  the  shock  front  (ignoring  any  low  level 

precursor5),  where  f(h)  is  a  fitted  function  solely  of  the 


3  The  same  procedure  applies  for  measured  radial  stress. 

4  Here  we  explicitly  distinguish  between  a  gage,  u(hQ,t), 
and  a  flow  surface,  u(h,t). 

5  Elastic  precursors  occurring  in  the  field  at  low 
stresses  introduce  the  problem  of  being  unable  to 
precisely  identify  f(h)  if  nothing  else  is  known.  The 
problem  is  exacerbated  if  the  signal  is  near  the  noise 
floor  and/or  discrete  resolution  of  the 
instrumentation.  This  imprecision  turns  out  to  be  a 
potentially  important  source  of  error  in  the  low  stress 
regime . 
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Lagrangian  gage  position;  suitable  functions  will  be  discussed 
later . 


Contrary  to  a  large  body  of  thought,  a  complete  stress  field 
and  motion  field  and  their  approximation  is  not  necessary  for 
estimation  of  the  strain  path  at  a  point.6  Indeed,  in  complex 
one-dimensional  flow,  seeking  the  average  stress  and  motion  flow 
across  a  large  expanse  of  the  field,  as  measured  in  several 
widely  spaced  locations,  is  an  improper  approach  since  strain  is 
a  quantity  defined  at  a  point.  In  practice,  what  can  actually  be 
calculated  are  relative  displacements  between  discrete  gages. 

One  would  ideally  like  to  have  gages  as  close  as  possible  (within 
limits)  to  the  point  at  which  the  path  is  sought  in  order  to  have 
the  average  calculated  strains  approach  the  actual  point  strains . 
This  is  especially  true  since  in  geological  media  it  is  nearly 
impossible  to  establish  the  regional  homogeneity  of  a  material 
undergoing  sharp  load-unloading. 

We  consider  the  following  definition  of  the  new  independent 
variables  x  and  y 

x  =  h; 

y=t-ta=t-  f ( h)  (5) 

and  define  the  transformation  u(h,t)  -*■  fl(x,y);  the  Jacobian, 

I J I  =  -1,  and,  thus,  if  H(x,y)  is  known  or  approximated,  we  are 
certain  we  can  always  obtain  the  global  inverse  transformation 
u  (  h ,  t )  «■  0  ( x,y )  . 

The  transformation  u(h,t)  -*  ft(x,y)  has  the  practical  effect 
of  aligning  the  characteristic  lines  of  (generally)  inhomogeneous 


It  is,  however,  a  necessity  for  calculation  of  stress 
bounds  and  is  helpful  in  gage  validation  efforts. 
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flow  so  that  they  all  begin  at  (x,0)  (see  Fig.  lb).  Exterior 
derivatives  across  the  lines  of  flow  in  {x,y}  can  then  be  much 
more  accurately  approximated  than  is  possible  in  {h,t}. 

We  note  the  following  identity 

u(h/t)  =  Ct[h,  t  -  f  (h)  ]  .  (6) 

The  transformation  (Eqn.  6)  introduces  no  inherent  error. 

However,  it  does  presume  that  f(h)  is  accurately  known.  This  is 
equivalent  to  the  presumption  that  the  initial  gage  locations  are 
accurately  known  and  that  the  gage  timing  is  correct. 


SECTION  4 

APPROXIMATING  THE  TRANSFORMED  FLOW  SURFACE  AT  X,Y 

The  transformed  surface  may  be  approximated  at  a  point  x0,y0 
as 

“(x.yj'x  v  =  “1(x'y0)  +  ^2{xo'y)  '  u(xo,yo)  (7) 

o' y  o 

and  globally  as  piecewise  continuous  polynomials.  The 
approximation  is  thus  a  numerical  surface  described  by  orthogonal 
splines . 

In  particular,  we  choose  Uj  and  u2,  quadratic  in  y  and 
quadratic  in  x,  as 

0l(xfyo)  =  ao  +  alx  +  a2*2'  (8) 

a2(xo'Y)  =  bo  +  biY  +  b2Y2-  (9) 

We  also  wish  to  consider  Q x  quadratic  in  ln(x)  or 

ln(a1(x,y0))  =  a0  +  axln(x)  +  a2ln2(x)  (10) 

in  order  to  more  accurately  approximate  a  flow  surface  that 
attenuates  very  steeply  with  increasing  range. 

The  function  u1(x,y0)  is  approximated  by  the  observed 
transformed  values  at  three  gages  at  a  constant,  transformed 
time,  yG.  Generally,  the  number  of  gages  available  at  distinct 
ranges  is  limited;  we  assume  here  that  three  gage  records  are 
available,  an  ’’inner’'  and  "outer"  gage  and  a  center  gage  at  the 
range  of  interest.  Thus,  the  approximating  quadratic  is 


perfectly  determined  and  is  given  by  an  augmented  Vandermonde 
determinant 
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(11) 


with  a  similar  expression  holding  for  ln(u1).  An  example  of  such 
a  fit  is  shown  in  Figure  2. 

The  situation  of  finding  the  three  constants,  bD,  bx,  and 
b2 ,  which  are  used  to  approximate  the  velocity  along  the  flow 
must  be  approached  by  a  different  route.  The  accuracy  of  locally 
approximated  derivatives  along  a  velocity  gage  can  be  severely 
affected  by  initial  sampling  rate,  resampling  rate,  noise  and 
record  band-width,  and  unexplainable  jerks  and  jumps.  (It  will 
turn  out  that  a  discontinuity  in  velocity  (leading  to  an  impulse 
in  acceleration)  will  be  integrated,  leading  to  a  unit-step 
discontinuity  in  calculated  strains  unless  smoothed  over.)  The 
result  is  that  while  a  low  order  least-squares  polynomial  spline 
is  used  to  locally  approximate  the  spline  u2,  the  fit  must  be 
made  carefully  in  the  sense  that  velocity  and  displacement 
recomputed  from  the  approximated  acceleration  record  must  be 
shown  to  be  consistent.  Both  linear  and  quadratic  splines  were 


tried . 

Figure  3a  shows  an  accelerometer  response  for  an  event  with 
its  integrated  velocity  record.  The  estimated  acceleration  using 
the  velocity  record  for  two  different  trials  is  shown  in  Figure 
3b.  The  sampling  rate  was  10  )isec  per  point,  which  insured 
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Figure  2. 


Orthogonal  splines  on  the  transformed  surface  used 
for  calculating  the  required  partial  derivatives  in 
the  neighborhood  of  the  point,  u(x0,y0).  The  splines 
are  approximated  in  different  ways. 
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Figure  3a.  Measured  acceleration  in  recompacted  dry  alluvium. 
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Figure  3b. 


Measured  acceleration  in  recompacted  dry  alluvium 
compared  to  two  estimates  from  velocity. 


VELOCITY  (m/s) 


adequate  resolution  of  the  peak.  Trial  1  used  a  linear  spline 
with  a  moving  window  of  20%  t ime-to-peak  until  peak  and  a  larger 
window  (by  a  factor  of  4)  during  decay.  Trial  2  used  a  quadratic 
spline,  again  with  a  moving  window  of  20%  time-to-peak;  the 
window  during  decay  was  multiplied  by  a  factor  of  10  in  this 
case . 

Both  fits  are  nearly  identical  except  for  the  error  in  Trial 
1  introduced  by  the  linear  spline  averaging  across  the  rather 
small  rise  in  velocity  immediately  after  peak.  The  quadratic 
spline  in  Trial  2  is  a  nearly  perfect  estimate  of  the  original 
accelerometer.  No  filtering  was  done  to  the  velocity  record  or 
the  estimated  accelerations  in  this  case.  Both  splines  appear  to 
be  adequate.  Experience  has  shown,  however,  that: 

1.  Equivalent  accuracy  with  a  linear  spline  requires 
smaller  windows  and  higher  sampling  rates  than  for  a 
quadratic  spline. 

2.  A  larger  window  during  decay  is  needed  when  using  a 
quadratic  spline  than  with  a  linear  spline;  this  is 
especially  true  with  decays  exhibiting  an  'oscillatory' 
character.  An  insufficiently  windowed  quadratic  spline 
has  the  nasty  tendency  to  capture  regions  of  inflection. 
This  leads  to  extremely  high  and  unrealistic  impulses  in 
estimated  acceleration. 

3.  Extremely  noisy  velocity  records  mandate  a  linear  spline 
(or  some  degree  of  filtering,  see  below). 

Ideally,  some  type  of  smoothing  should  be  done  to  allow 
small  fitting  windows  to  be  used  so  that  the  detailed  structure 
of  the  velocity  waveform  is  captured  (not  attributable  to 
stochastic  processes).  A  finite  impulse  response  (FIR)  low  pass 
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filter  with  linear  phase  characteristics  that  is  useful  is  an  n- 
repeated  three  point  moving  average  filter  (Ref.  13); 


u  . 
1 


+  u  . 
] 


(12 


The  filter  usually  only  needs  to  be  applied  after  rise-to-peak; 7 
the  number  of  applications  to  a  record  depends  on  the  initial 
band  width  and  power  spectra.  This  filter  is  especially  useful 
for  effectively  filtering  calculational  records,  which  generally 
s[Aw  clear  trends  despite  the  presence  of  numerical  oscillations 
about  the  trend. 


The  untransformed  derivatives  for  the  approximating  splines 
u i  and  u2  given  by  Equations  8  and  10  are  derived  next.  The 
starting  point  is  the  earlier  recognition  that  u(h,t)  =  u[h,t  - 
f(h)j.  Thus,  the  spatial  and  temporal  derivatives  are,  in  terms 
of  x  and  y. 
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where  t  -  f(h)  is  being  treated  as  a  new  independent  variable, 
say  z.  Using  the  chain  rule. 


The  reason  being  that  early-time  accelerations  are  very 
high  and  easily  approximated  with  sufficient  sampling 
of  the  velocity  waveform .  Accelerations  after  peak 
velocity  can  easily  be  several  orders  of  magnitude 
smaller  than  the  initial  accelerations .  (See  Fig.  3a.) 
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In  problems  like  this  it  is  often  helpful  to  explicitly  indicate 
which  independent  variable(s)  is  (are)  being  treated  as  constant, 


The  approximating  surface  in  {x,y}  is  given  by 


ao  a1+a2ln(x)  2  „ 

u  ( x ,  y )  =  e  x  +  b0  +  +  *>2y  -  u(xQ,yo).  (16) 


Therefore, 


u(h,t)  =  u ( x, z )  = 

ao  al+a2ln(x) 


e  x 


+  b  +  b.z  +  b_ z  -  u ( x  , y  ) . 


(17) 


Hence, 


=  ea°[ai  +  2a2ln(x)]x 


a1~l+a2ln(x) 


,  and. 


(18) 


Q,u_Lx^j-l~|  =  h  +  2b0z . 


dz 


(19) 
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And  finally,  substituting  these  derivatives  into  Equation  15: 


8u(h*tn 
9h  Jt 


a  a.-l+a?ln(h) 

e  °[a1  +  2a2ln(h) ]h 


-  fh(h)[b1  +  2b2(t  -  f (h) ) ] . 


(20) 


To  evaluate 

the  acceleration 

in 

{ h , t }  we 

again 

use  Equation 

13 

5u  l  h ,  t ) 
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This  reduces  to 

— atzt))h  =  bi  +  2b2[t  "  f(h>l-  (22) 

And,  of  course,  acceleration  is  seen  to  be  invariant  whether  it 
is  approximated  in  {x,y}  or  {h,t}. 

Important  to  note  is  that  the  preceding  approximations  are 
not  global;  they  are  valid  only  in  the  region  {x0,y0}. 

Therefore,  while  the  approximating  function  may  be  legitimately 
differentiated  at  x0,y0,  integration  outside  the  immediate 
neighborhood  of  the  point  is  not  permissible. 

Nothing  has  been  said  of  the  form  of  f(h)  or  means  for  its 
approximation.  Three  different  functional  forms  found  useful  are 
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M 


ta  =  f(n)  = 


aG  +  axh 


=  a„  +  a-,h  +  a,h2 


ao  al+a2ln^ 
=  e  h 


Selection  and  fitting  of  a  form  of  f(h)  can  be  done  in  two  ways. 
In  one  approach,  various  forms  can  be  fitted  to  all  data  points. 
Then  the  form  which  minimizes  error  according  to  any  of  several 
criteria  is  used.  The  second  approach  is  to  fit  f(h)  locally, 
using  data  points  in  the  neighborhood  of  the  differencing  gage 


location . 


The  most  appropriate  approach  is  problematical.  On  the  one 
hand,  if  gage  time-of-arrivals  (or,  especially,  surveyed  gage 
positions)  are  highly  suspect,  the  only  possible  approach  would 
appear  to  be  to  estimate  f(h)  in  the  neighborhood  by  a  fit  to  all 
available  data  points.  On  the  other  hand,  if  gage  positions  and 
time-of-arrivals  are  known  to  be  accurate,  approximation  of  f(h) 
should  be  done  in  the  neighborhood  of  h0. 


SECTION  5 

APPROXIMATING  STRAIN-TIME  HISTORIES 


We  start  with  the  following  definitions  of  strain  in  one¬ 
dimensional  spherical  flow:8 
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(24) 


Note  that  ev  is  the  full  strain  of  a  unit  volume  element  with 
principal  axes,  ^r,  Xd  >  Xt>  undergoing  extensional  strains  of  -er 
and  -et  =  eg.  That  is. 


e  =1-  n  (1+e.) 
v  i  =  ! 


In  order  to  explicitly  evaluate  the  above  strains  in  terms  of  the 
approximate  surface,  we  begin  with  the  Eulerian  position 
expressed  in  terms  of  the  Lagrangian  gage  position  and  the  gage 
response  as  a  function  of  Lagrangian  position  and  time. 

Using  Equation  4,  we  can  write  the  strains  as: 


1 

€t(ho't)  =  ~h"~  *  J  u(h  ,t)dt, 

o  f(hJQ)  ° 

er(ho  t)  =  -  J  9U9h't)dt> 

r  o,  f(h')  on 


(25) 


A  strain  which  produces  compression  is  positive  in  this 
text . 
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•dt  . 


f(ho) 


Other  commonly  employed  variables,  p  and  p  are  found  from  ev  by 
the  following  formulas: 

e  p 

p  «  y~TT~  and  P  =  •  ( 2 6) 

v  v 

Generally,  the  number  of  long-term  surviving  gages  providing 
acceptable  velocity-time  histories  is  small  enough  that  a  global 
fit  to  the  entire  flow  surface  is  an  ill-defined  undertaking. 
Practically,  Equations  25  are  best  evaluated  in  discrete  time- 
steps  with  a  local  approximating  surface,  u(h,t)  found  for  each 
time-step . 

The  technique  resulting  from  this  paper  assumes  at  least 
three  gages9  are  available  and  uses  the  approximating  polynomials 
from  the  previous  section  to  calculate  the  spatial  and  time 
derivatives  at  each  time-step.  The  defining  equations  for 
strains  are  integrated  over  each  time-step  and  the  terms 
cumulatively  summed.  A  typical  term  using  a  quadratic  polynomial 
(Eq.  8)  fit  across  the  flow  {h  direction)  and  a  quadratic  fit 
(Eq.  9)  along  the  flow  (acceleration  as  a  function  of  time)  is: 


f(h  )+(k+l)At  „ 

°  r  auth.t  -  f ( h) ] 

f(hJo)+kAt  9h 


dt  =  Atfa^)  +  2h  ai 

o  2 


(M  _  h(k) 


Use  of  a  larger  number  of  gages  to  approximate  the  a^s 
is  certainly  possible  and,  in  some  cases  (i.e.,  if  the 
strain  field  is  known  a  priori  to  be  exceptionally 
uniform),  probably  preferable.  However,  as  mentioned 
before,  there  are  definite  reasons  why  a  field-wide  fit 
should  not  be  done. 


•  fh(ho)]  -  A2t  •  (2k  +  1)  •  fh(hQ)  •  b<k),  (27) 

where  each  of  the  and  the  is  approximated  at  f(hQ)  +  kAt  as 
indicated  previously.  The  second  order  term  is  a  result  of  the 
quadratic  fit  to  the  velocity  and  cannot  be  neglected.  Also, 
note  that  the  limits  of  integration  are  evaluated  at  a  specific 
gage,  hQ,  and  are  considered  constant  for  each  integration. 

Thus,  for  example,  a  numerical  solution  to  Equation  25  (er) 
is  derived  in  full 
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If  Equation  10  is  used  to  approximate  the  flow  across  the  gages 
instead  of  8,  the  solution  is 
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A2t  •  (2k  +  1)  •  fh(hQ)  •  b^k; 


where  kmax  corresponds  to  the  index  of  the  final  time-point.  The 


ratio,  |t-J  in  Equation  24  is  equal  to: 


(30) 


where  the  integral  must  be  evaluated  numerically.  Note  that  the 
integral  is  simply  a  displacement  history  along  a  gage  (h 
constant)  and  may  be  evaluated  by  several  well  known  schemes. 
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SECTION  6 

CALCULATION  OF  STRESS  DIFFERENCE  FROM 
COMBINED  VELOCITY/STRESS  RECORDS 


The  conservation  of  momentum  equation  can  be  solved  for  the 
stress  difference  as  follows: 
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The  spatial  derivative  of  radial  stress  can  be  approximated  in 
the  same  manner  as  was  done  for  velocity;  the  indicated  strains 
are  presumably  already  at  hand.  Thus,  with  three  velocity 
records  and  three  stress  records,  the  stress  difference  at  the 
location  of  the  center  gage  can  be  found.  The  location  of  the 
center  stress  and  velocity  gage  must  coincide.  If  they  do  not, 
interpolation  across  the  gages10  can  be  done  to  provide  a  "gage." 

Calculation  of  x  hinges  on  accurate  determination  of 
fi(hQ,t),  and  use  of  the  splines  and  filtering  methods  of  Section 
4  will  do  a  lot  to  reduce  a  noisy  result  to  a  clear  stress  path. 
In  practice,  estimation  of  ft  will  be  the  limiting  factor  in 
estimating  x. 

Note  that,  to  first  order,  the  stress  difference  is  in  error 
directly  proportional  to  the  error  in  displacement  due  to  the 
term ,  r ( h , t ) . 
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See  Section  4 . 


SECTION  7 
STRESS  BOUNDS 


Trulio's  (Ref.  6)  lower  stress  bound  equation  can  be 
obtained  directly  from  the  equation  of  momentum  conservation  by 
assuming  the  media  is  a  perfect  fluid  and,  hence,  the  stress 
difference  is  identically  equal  to  zero.  Then  integrating  along 
a  path  of  constant,  t, 
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(32) 


Figure  4  is  given  as  an  aid  in  visualizing  the  paths  of 
integration  in  this  integral  and  the  one  that  follows. 

The  corresponding  Lagrangian  integral  is  found  by  observing 
that  the  total  differential,  dr,  is 

dr  -  If) dh  +  If] dt  <33» 


and  along  a  path  of  constant  time: 
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(34) 


where  the  final  integral  in  Lagrangian  coordinates  is  found  by 
applying  conservation  of  mass  ( Eq .  1)  to  eliminate  9r/9h.  Note 
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Figure  4. 


Surface  of  lower  stress  bound  integrand  using  spherical  cavity  test 
problem.  The  lower  stress  bound  integral  is  given  by  the  area  under 
the  displayed  surface  along  a  path  of  constant  time  projected  onto  the 
■x,y  plane,  beginning  at  h0  and  ending  at  hf.  hf  is  a  function  of 
time,  making  the  stress  bound  a  function  of  h0  and  time.  Below,  the 
path  of  constant  time  (a  linear  function  in  the  test  problem)  is  shown, 
along  with  the  gage  line,  h0,  at  a  range  of  1.5  and  time  of  7.3,  and 
the  upper  limit  of  the  integral  at  hf.  Note  that  hf  is  at  the  advancing 
wave  front  at  y  =  0;  further  integration  along  the  path  would  contribute 
nothing  to  the  integral.  Note  that  integration  at  ranges  where  hf 
exceeds  the  measured  field  (h  -  10,  here),  the  stress  bound  requires 
extrapolation.  The  possible  error  is  not  always  small. 


that  the  upper  limits  of  integration  are  identical  since  r  =  h  at 
the  advancing  wave  front. 


Under  the  assumption  that  tangential  stresses  can  never 
exceed  the  material  tensile  limit,  8,  an  upper  bound  to  radial 
stress,  can  be  found,  viz.: 
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Evaluation  of  Equations  34  and  35,  which  are  line  integrals 
over  the  transformed  surface  (that  is,  the  path  of  integration  is 
not  along  either  of  the  transformed  axis),  is  somewhat  subtler 
than  the  estimation  of  strains,  which  merely  requires 
approximation  of  the  local  function  surface  at  a  point. 

Some  difficulty  is  also  introduced  in  the  following  way: 
the  functional  surfaces  used  in  the  strain  relations  are  those  of 
velocity  and  displacement,  surfaces  which  tend  to  be  smooth 
(especially  the  displacement  surface),  and  numerically 
predictable . 

The  two  surfaces  described  by  the  integrands  of  the  stress 
bound  formula  are  rather  complex,  being  functions  of  both 
acceleration  and  reciprocal  displacement.  We  have  found  the 
following  numerical  technique  to  give  answers  which  match  well 
with  known  results. 

The  profiles  of  the  transformed  surface  to  be  approximated 
are  calculated  per  the  methods  of  Section  II.  A  global 
approximation  to  the  surface  is  then  constructed  using  the 
following  basis  functions  (Refs.  14,  15, and  16): 
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SECTION  8 

TEST  PROBLEMS  FOR  SOLUTION  METHODS 


Several  different  numerical  solutions  of  the  equations  of 
spherical  flow  in  soils  and  rocks  have  been  reported  (Refs.  1,  3, 
and  4).  Generally,  the  methods  have  first  been  applied  to 
waveforms  from  problems  with  known  solutions  in  order  to  obtain 
an  idea  of  the  errors  and  difficulties  associated  with  each 
numerical  method.  This  is  done  here. 

The  following  analytic  problem  is  adapted  in  large  part  from 
Grady,  (Ref.  1),  and  was  used  to  ascertain  the  accuracy  of  our 
solution  procedure.  In  addition,  the  method  is  applied  to 
ascertain  the  accuracy  of  the  solution  method  used  in  computing 
the  stress  bounds  introduced  by  Trulio.  In  addition,  stress  and 
velocity  waveforms  from  a  1-D  spherical  calculation,  simulating  a 
contained  nuclear  detonation  (yield  =  15.15  kt)  in  Yuma  alluvium 
were  used  as  more  representative  of  the  input  for  which  the 
method  was  developed. 

Consider  an  expanding  cavity  (with  prescribed  elastic 
boundary  conditions)  in  an  infinite  homogeneous  isotropic  medium. 
The  field  solution  can  be  written  in  terms  of  a  scalar  elastic 
potential;  for  waves  moving  away  from  the  cavity  the  solution  in 
terms  of  dimensionless  length  (H)  and  time  (T)  is 


(36) 


where  U  is  the  displacement  as  seen  by  a  fixed  observer  in  an 
inertial  (Eulerian)  frame.  We  choose  v  =  0.25,  K  =  5/9  and  G  = 
1/3,  and  the  dimensionless  rest  density  and  elastic  sound  speed 
will  equal  unity.  The  boundary  conditons  are 


Displacement  at  cavity:  U(1,T)  =  g(l  -  T),  T  >  1, 
Displacement  vanishes  at  wave  front:  U(T,T)  =  0,  T  >  1. 

Let  £  =  H  -  T.  The  solution  in  terms  of  the  cavity 
displacement  is: 


f(S) 


e^  g(£' )d$' 


(37) 


All  field  quantities  of  interest  can  be  calculated  in  terms  of 

f  (£)  / 


U  =  -H~2  f  +  H  1 f 


V  =  H~2  f  '  -  H~ 1 f 


=  -H-2  f '  '  +  H~ 1 f '  '  ' 


ev  -  -H-h' 


er  =  -H  1f"  +  2H-2 f '  -  2H-3 f 
et  =  -H-2  f '  +  H-3  f 
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+  (1 


(38) 


If  we  describe  the  same  test  problem  as  Grady,  g(£)  =  0 . 02£2e£, £ 
<  0 ,  then 


f  =  0.00667£3e£, 

f '  =  £2e£(0.02  +  0.00667£), 

f"  =  £e£(0. 00667^2  +  0.04^  +  0.04),  (  39) 

and  f'"  =  0.00667^2  +  0.06£  +  0.12  +  0.04s-1). 

Figure  5  is  the  results  of  the  method  applied  to  the  problem  of 
estimating  the  strains  at  H  =  6 ,  using  the  known  velocity 
waveforms  at  H  =  5 ,  6,  7  (Eq.  38).  At  a  sufficiently  high 
sampling  rate,  the  difference  between  the  known  strains  and  the 
estimated  strains  is  very  small. 

The  upper  and  lower  stress  bounds,  in  terms  of  the  test 
problem,  are: 
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where  the  final  integral  for  the  lower  bound  is  an  approximation 
derived  by  assuming  U  <<  H.  This  is  generally  a  valid  assumption 
at  moderate  to  distant  ranges  and  early  to  moderate  times.  The 
integrals  are  evaluated  numerically. 


Results  of  the  method  for  estimating  stress  bounds  discussed 
in  Section  VI,  compared  to  known  upper  and  lower  solutions 
calculated  at  H  =  5  from  Equations  40  and  41,  are  shown  in  Figure 


Figure  5.  Comparison  of  Lagrangian  strain  estimation  with  analytical  spherical  cavity 
test  problem. 


6a.  Surface  profiles  were  calculated  at  H  =  6,  7,  8.  The 
tensile  limit  was  zero  in  Equation  41.  The  late-time  match  is 
very  good.  This  is  to  be  expected,  as  the  length  of  the 
integration  interval  grows  and  encompasses  more  profiles.  Error 
at  early  time  is  to  be  expected  for  two  reasons: 

1.  The  known  accelerations  in  the  test  problem  have  a 
finite  value  at  t0+. 

2.  The  step  size  is  necessarily  finite. 

However,  the  errors  are  quite  small. 

Figure  6b  shows  the  magnitude  of  the  integrands  at  three 
times--early ,  moderate,  and  late,  as  a  function  of  H.  This  plot 
provides  some  indication  of  where  the  largest  errors  might  be 
expected  to  occur. 

Note  that,  at  late  times,  increasingly  large  contributions 
arise  from  the  accelerations  occurring  at  the  surface  close  to 
Hf.  In  practical  terms,  the  late-time  stress  bounds  must  include 
an  error  due  to  any  indeterminancy  in  measuring  velocity  rise 
times.  That  is,  rise  times  that  are  too  long  will  give  a  falsely 
lower  bound. 

A  one-dimensional  finite  difference  calculation  was  run  in 
order  to  provide  a  test  loading  more  comparable  (in  both  stress 
and  material)  to  loadings  of  current  interest. 

Pertinent  information  about  the  calculations: 
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Figure  6a.  Comparison  of  estimated  stress  bounds  with  analytical  s 
cavity  test  problem. 
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5.5  6.0 

DIMENSIONLESS  RANGE 


4.  Point  of  interest  (center  point):  20  meters 

5.  Peak  velocity  at  POI :  Approximately  650  m/sec 

6.  Peak  radial  stress  at  POI:  Approximately  30  kbar 


Figure  7  shows  the  stress  and  velocity-time  histories  from 
the  calculation. 

Figure  8  shows  the  estimated  radial  and  volumetric  strains 
compared  against  the  known  strains  from  the  calculation. 

The  estimated  strains  initially  appear  to  match  well  with 
the  known  strains  though  with  decreasing  accuracy  as  the  width  of 
the  estimation  interval  increases. 

Considering  Trial  3  to  be  grossly  in  error,  we  would  predict 
that,  to  first  order,  estimation  intervals  in  this  media  at  these 
initial  loading  conditions  should  be  limited  to  the  interval  in 
Tr ia 1  2 : 

— — - — — — =  6.0  m/kt1^  (fully  contained) 

15.15  kt 

The  effect  of  the  error  in  Trial  1,  while  seemingly  small  at 
the  scale  shown,  is  actually  enormous,  when  considered  at  the 
scale  that  would  be  required  for  estimation  of  the  material 
model . 

Figure  9  shows  the  known  strain  paths  in  the  calculation  at 
20  meters  compared  with  the  strain  path  obtained  from  the 
Lagrangian  analysis,  the  stress  used  being  the  calculated  stress 
and,  thus,  not  subject  to  validation.  The  estimated  model  is 
clearly  diametrically  opposite  in  all  respects  to  the  model  that 
was  actually  used. 


for  differences  in  loading. 


The  source  of  error  is  due  entirely  to  the  discrete  zone 
size  used.  This  allowed  velocity  (and,  hence,  strains  estimated 
from  velocity)  to  rise  to  peak  more  quickly  than  compression, 
which  required  the  shock  front  to  traverse  the  zone  completely. 
In  this  case,  the  difference  in  timing  that  produced  such  a 
completely  invalid  material  model  estimation  was  60  /i sec,  or 
about  ten  calculational  cycles. 

The  lesson  learned  above,  that  small  errors  in  estimated 
strains  lead  to  extremely  large  errors  in  model  estimation, 
applies  in  practice  in  a  somewhat  different  manner.  We  can  note 
that  gage  rise  times  play  a  singular  role  in  determining  the 
accuracy  of  an  estimated  material  model. 
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SECTION  9 


l.C»  ■  4 


SENSITIVITY  OF  ev  TO  BASELINE  CORRECTIONS 


APPLIED  TO  VELOCITY  RECORDS 


It  is  important  to  investigate  the  effect  of  baseline 


corrections  because: 


1.  Some  sort  of  velocity  drift  correction  is  inevitably 
applied  to  integrated  accelerometer  records,  and 


2.  Such  a  correction  implies  a  direct  effect  on  so-called 
observed  radial  displacements ,  and  et,  er,  and  ev  are 
all  functions  of  this  displacement  and  partials  thereof 


Consider  ev  at  some  fixed  time  after  arrival  of  motion. 


Then,  by  definition: 


e  =  1  -  f~1 

v  LhJ 


9r(h, t 
9h 


Then  write  r  as: 


r ( h , t  )  =  h  +  u(h,t)dt  =  h  +  D(h,t  ), 

^f(h)  ° 


where  D  is  the  displacement  at  the  fixed  time, 
function  of  the  independent  variables  h  and  D. 


Thus,  ev  is  now  a 


Then : 


!  ,  x  3D  ,  D 

e  =1-  1  +  -=rr  •  1  +  u 

v  9hJ  l  hj 


Changes  in  ev  due  to  small  changes  in  h  (due  to  positional 
errors)  and  D  (due  to  baseline  shifts)  are  given  by 


4ev  =  40  *  §D(£v(h,Dnh  D  +  4h  •  fh<£v(h'D)>h  ,D  •  <3 

o  o  o  o 


This  is  the  equation  of  the  tangent  plane  at  (h0,D0). 
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(Ah  =  0) 


This  expression  is  difficult  to  evaluate  fully.  In  order  to 
evaluate  it  qualitatively,  assume  1  +  D0/h0  is  approximately  1, 
which  it  will  be  for  early  times.11  Then: 


(4 


11  This  will  prove  to  be  a  valid  operation,  since  these 
terms  are  greater  than  zero.  It  also  removes  the 
stipulation  for  earlier  times. 


40 


£ & 


Therefore , 1 2 


ev  -  evo  =  (D  -  Dc)  •  (positive  number) 

D  >  D0  -*•  ev  >  e  vo ,  (42) 

that  is,  a  baseline  shift  that  gives  higher  displacements  for  the 
central  gage  at  some  time  works  to  give  a  greater  volumetric 
strain  at  the  same  time.  An  example  of  this  type  of  error  is 
shown  in  Figure  10.  Other  baseline  shifts  for  the  two  end  gages 
can  also  be  calculated. 

This  result  has  the  important  implications  that  positive 
baseline  correction  (a  correction  which  gives  a  high 
displacement)  works  to  keep  ev  (and  hence  /x)  higher  in  magnitude 
than  it  actually  is.  Then,  the  strain-time  history,  in 
conjunction  with  an  accepted  stress-time  record,13  would  tend  to 
indicate  an  abnormally  high  unload  modulus  or  even  a  negative 
unload  modulus,  "competency,"  or  the  continued  compaction  of 
material  with  falling  pressure. 

A  drift  correction  of  opposite  magnitude  (giving 
displacements  that  are  too  small)  would  tend  to  indicate 


The  term  92h/9D2  is  asymptotic  to  and  always  greater 
than  zero;  at  early  times  it  will  be  very  large.  This 
can  be  seen  by  noting  that  the  graph  of  h  versus  D  in  a 
divergent  field  will  be  a  power  law  curve  upwardly 
concave.  9h/9D  will,  therefore,  always  be  negative  and 
upwardly  convex.  Note  also  that  3D/3h  will  be  a 
relatively  large  negative  number;  hence,  the  first  term 
in  the  last  parentheses  is  neglected.  The  statement 
follows  easily  from  these  arguments. 

An  assumption  here  being  that  the  times-of -arrival  of 
the  stress  and  velocity  records  are  identical. 
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excessive  "dilatancy"  during  stress  relief--too  much  unload 
recovery . 

The  following  two  statements  are  especially  pertinent  in 
materials  that  unload  very  little,  such  as  dry  soils,  where  the 
slope  of  volumetric  strain  is  very  shallow  with  respect  to  time. 

Linear  baseline  corrections  are  inevitably  applied  upon 
shock  passage;  if  incorrect,14  such  corrections  thus  give 
resultant  displacements  that  become  increasingly  in  error  with 
the  passage  of  time.  Volumetric  material  behavior  on  unload 
derived  from  Equation  24  is,  thus,  immediately  subject  to  a  doubt 
that  grows  with  the  degree  of  unloading.  The  degree  of  error  is 
unknown  unless  the  velocity  records  are  known  a  priori  to  be 
correct.  Tentatively  correct  corrections,  of  course,  imply 
tentatively  correct  strains  which,  therefore,  must  be  verified  by 
other  methods. 

It  might  also  be  noted  that  "averaging"  records  in  a 
neighborhood  in  an  attempt  to  cull  a  more  accurate  record  is 
incorrect  unless  it  is  known  that  the  distribution  of  gage  drifts 
is  uniformly  random  about  zero.  (For  example,  an  average  of 
three  records  in  a  neighborhood  cannot  produce  a  better 
approximation  than  any  one  of  the  individual  records  if  all  three 
gages  require  systematic  positive  drift  corrections  during 
unload . ) 


14 


Frequently,  more  than  one  correction  is  possible, 
especially  when  like-gage  responses  are  compared. 
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ATTN:  ATSB-CTD 
ATTN:  TECH  LIBRARY 

U  S  ARMY  BALLISTIC  RESEARCH  LAB 
ATTN:  SLCBR-DD-T 
ATTN:  SLC8R-SS-T  (TECH  LIB) 

U  S  ARMY  COMB  ARMS  COMBAT  DEV  ACTY 
ATTN:  ATZL-CAP 

U  S  ARMY  COMD  &  GENERAL  STAFF  COLLEGE 
ATTN:  ACQ  LIBRARY  DIV 
ATTN  ATSW-TA  D 

U  S  ARMY  EUROPE  AND  SEVENTH  ARMY 
ATTN:  AEAGC 

U  S  ARMY  INFANTRY  CTR  &  SCH 
ATTN  ATSH  CD  CS 

U  S  ARMY  MATERIEL  COMMAND 
ATTN:  AMCCN 

U  S  ARMY  MATERIEL  SYS  ANALYSIS  ACTVY 
ATTN:  AMXSY-CR 

U  S  ARMY  M'SSILE  COMMAND 
ATTN:  AMSMI  XF 

U  S  ARMY  NUCLEAR  &  CHEMICAL  AGENCY 
ATTN:  LIBRARY 

U  S  ARMY  STRATEGIC  DEFENSE  CMD 
ATTN:  DACS  BM/J  KAHLAS 

U  S  ARMY  TRADOC  SYS  ANALYSIS  ACTVY 
ATTN:  ATRC  TO 

U  S  ARMY  WAR  COLLEGE 
ATTN:  LIBRARY 

USA  MILITARY  ACADEMY 

ATTN  DOCUMENT  LIBRARY 

V  CORPS 

ATTN  COMMANDER 
ATTN  G  3 

VII  CORPS 

ATTN  G  3 


DEPARTMENT  OF  THE  NAVY 

CARRIER  AIRBORNE  EARLY  WARNING  WING  12 
ATTN:  COMMANDER 

CARRIER  GROUP  1 

ATTN:  COMMANDER 

CARRIER  GROUP  2 

ATTN  COMMANDER 

CARRIER  GROUP  3 

ATTN:  COMMANDER 

CARRIER  GROUP  4 

ATTN:  COMMANDER 

CARRIER  GROUP  5 

ATTN  COMMANDER 

CARRIER  GROUP  7 

ATTN:  COMMANDER 

CARRIER  GROUP  8 

ATTN:  COMMANDER 

CRUISER  DESTROYER  GROUP  ONE 
ATTN:  COMMANDER 

CRUISER  DESTROYER  GROUP  12 
ATTN  COMMANDER 

CRUISER-DESTROYER  GROUP  2 
ATTN:  COMMANDER 

CRUISER  DESTROYER  GROUP  3 
ATTN:  COMMANDER 

CRUISER-DESTROYER  GROUP  5 
ATTN:  COMMANDER 

CRUISER  DESTROYER  GROUP  8 
ATTN:  COMMANDER 

DAVID  TAYLOR  NAVAL  SHIP  R  &  D  CTR 
ATTN:  CODE  174 

FIGHTER  AIRBORNE  EARLV  WARNING  WING 
ATTN:  COMMANDER 

FIGHTER  WING  1 

ATTN.  COMMANDER 

FLEET  INTELLIGENCE  CTR  EUROPE  &  ATLANTIC 
ATTN:  LIBRARY  CODE  113 

LIGHT  ATTACK  WING 

ATTN:  COMMANDER 

MARINE  CORPS  DEV  &  EDUCATION  COMMAND 
ATTN:  COMMANDER 

MEDIUM  ATTACK  TACTICAL  ELECTRONIC 
ATTN:  COMMANDER 

NAVAL  AIR  FORCE 

ATTN:  COMMANDER 
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NAVAL  AIR  FORCE 

ATTN:  COMMANDER 

NAVAL  INTELLIGENCE  SUPPORT  CTR 
ATTN:  NISC-30 

NAVAL  OCEAN  SYSTEMS  CENTER 

ATTN:  CODE  9642/TECH  LIB 

NAVAL  POSTGRADUATE  SCHOOL 

ATTN:  CODE  1424  LIBRARY 

NAVAL  RESEARCH  LABORATORY 
ATTN:  CODE  1220 
ATTN:  CODE  2627/TECH  LIB 

NAVAL  SEA  SYSTEMS  COMMAND 
ATTN:  SEA-09G53  (LIB) 

ATTN:  SEA-643 

NAVAL  SURFACE  FORCE 

ATTN:  COMMANDER 

NAVAL  SURFACE  FORCE 

ATTN:  COMMANDER 

NAVAL  WEAPONS  EVALUATION  FACILITY 
ATTN:  CLASSIFIED  LIBRARY 

NUCLEAR  WEAPONS  TNG  GROUP,  ATLANTIC 
ATTN:  CODE  221 
ATTN:  DOCUMENT  CONTROL 

NUCLEAR  WEAPONS  TNG  GROUP,  PACIFIC 
ATTN:  DOCUMENT  CONTROL 

OFC  OF  THE  DEPUTY  CHIEF  OF  NAVAL  OPS 

ATTN:  OP  654/TRAT  EVAL  &  ANAL  BR 
ATTN:  OP  955(AAW  DIV) 

ATTN:  OP  981 

OFFICE  OF  THE  CHIEF  OF  NAVAL  OPERATIONS 
ATTN:  CNO  EXEC  PANEL/OP  00K 

OPERATIONAL  TEST  &  EVALUATION  FORCE 
ATTN:  CODE  80 

OPERATIONAL  TEST  &  EVALUATION  FORCE, 
ATTN:  INTEL  OFFICER 

SUBMARINE  FORCE 

ATTN:  COMMANDER 

SUBMARINE  GROUP  6 

ATTN,  COMMANDER 

TACTICAL  TRAINING  GROUP,  PACIFIC 
ATTN  COMMANDER 

TACTICAL  WINGS  ATLANTIC 
ATTN  COMMANDER 

THEATER  NUCLEAR  WARFARE  PROGRAM  OFC 
ATTN  PMS423 

U  S  NAVAL  FORCES,  EUROPE 

ATTN  N54,  NUC  WARFARE  OFCR 


U  S  NAVY  SEVENTH  FLEET 

ATTN:  COMMANDER 

U  S  NAVY  SIXTH  FLEET 

ATTN:  COMMANDER 

DEPARTMENT  OF  THE  AIR  FORCE 

AF/INE 

ATTN:  INA 

AFIS/INT 

ATTN:  INT 

AIR  FORCE  SYSTEMS  COMMAND 
ATTN:  DL 
ATTN:  SD 

AIR  FORCE  WEAPONS  LABORATORY 
ATTN:  SUL 

AIR  UNIVERSITY  LIBRARY 
ATTN:  AUL-LSE 

ASSISTANT  CHIEF  OF  THE  AIR  FORCE 
ATTN:  SAF/ALR 

DEPUTY  CHIEF  OF  STAFF/AF  RDQM 
ATTN:  AF/RDQI 

DEPUTY  CHIEF  OF  STAFF/XOO 
ATTN:  AF/XOC 

DEPUTY  CHIEF  OF  STAFF/XOX 

ATTN:  AFXOXFM/PLNS,  FRC  DEV  MUN  PLNS 
ATTN:  AFXOXFS/FRC  DEV,  STRAT  OFF  FRC 

SPACE  DIVISION/CWH 

ATTN:  CWH  (DSCS  111) 

STRATEGIC  AIR  COMMAND/NRI  STINFO 
ATTN:  NRI/STINFO 

STRATEGIC  AIR  COMMAND/STIC 
ATTN:  STIC  (544SIW) 

TACTICAL  AIR  COMMAND  TAC/DOA 
ATTN:  TAC/DOA 

TACTICAL  AIR  COMMAND/XPJ 
ATTN:  TAC/XPJ 

U  S  AIR  FORCES  IN  EUROPE/DEX 
ATTN:  USAFE/DEXX 

U  S  AIR  FORCES  IN  EUROPE/DOT 
ATTN:  USAFE/DOQ 

U  S  AIR  FORCES  IN  EUROPE/INAT 
ATTN  USAFE/INAT 

2ND  ACCS 

ATTN  DOC 
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DEPARTMENT  OF  ENERGY 

DEPARTMENT  OF  ENERGY 
ATTN  OMA.  DP  22 

LAWRENCE  LIVERMORE  NATIONAL  LAB 
ATTN  L  35 
ATTN.  L  38 
ATTN  L-389 
ATTN  L  450  W  HOGAN 
ATTN  L  53  TECH  INFO  DEPT  LIB 
ATTN  Z  DIVISION  LIBRARY 

lOS  ALAMOS  NATIONAL  LABORATORY 
ATTN  D  RICHMOND 
ATTN  F601  T  DOWLER 
ATTN  MS  P364  REPORT  LIBRARY 
ATTN  R  SANDOVAL 

SANDIA  NATIONAL  LABORATORIES 

ATTN  TECH  LIB  3141/RPTS  RCVGCLRK 
ATTN  0333  R  B  STRATTON 

OTHER  GOVERNMENT 

CENTRAL  INTELLIGENCE  AGENCY 

ATTN  N10  STRATEGIC  SYS 

FEDERAL  EMERGENCY  MANAGEMENT  AGENCY 
ATTN.  GORRELLNPCP 
ATTN.  OFC  OF  RSCH/NP  H  TOVEY 

U  S  DEPARTMENT  OF  STATE 
ATTN  PM/TMP 

DEPARTMENT  OF  DEFENSE  CONTRACTORS 

ADVANCED  RESEARCH  &  APPLICATIONS  CORP 
ATTN  DOCUMENT  CONTROL 

BDM CORP 

ATTN.  CWASAFF 
ATTN:  J  80DE 
ATTN.  J  BRADDOCK 
ATTN:  R  BUCHANAN 

BOEING  CO 

ATTN.  H  WICKLEIN 
ATTN:  J  W  RUSSELL 

CALIFORNIA  RESEARCH  AND  TECHNOLOGY,  INC 

2  CYS  ATTN:  E  RINEHART 

2  CYS  ATTN  G LLOYD 

COMPUTER  SCIENCES  CORP 
ATTN:  FEISENBARTH 

GRUMMAN  CTEC,  INC 

ATTN  SSHRIER 

HORIZONS  TECHNOLOGY,  INC 
ATTN.  J  PALMER 

I  IT  RESEARCH  INSTITUTE 

ATTN  DOCUMENTS  LIBRARY 


INSTITUTE  FOR  DEFENSE  ANALYSES 
ATTN:  CLASSIFIED  LIBRARY 
ATTN:  JGROTE 

JAYCOR 

ATTN:  R  SULLIVAN 

KAMAN  SCIENCES  CORP 
ATTN:  F SHELTON 

KAMAN  SCIENCES  CORP 
ATTN:  E CONRAD 

KAMAN  SCIENCES  CORPORATION 
ATTN:  DASIAC 

KAMAN  TEMPO 

ATTN:  DASIAC 

MARTIN  MARIETTA  DENVER  AEROSPACE 
ATTN:  J  DONATHAN 

ORION  RESEARCH  INC 

ATTN:  J  E  SCHOLZ 

PACIFIC-SIERRA  RESEARCH  CORP 

ATTN:  H  BRODE,  CHAIRMAN  SAGE 

PACIFIC-SIERRA  RESEARCH  CORP 
ATTN:  D  GORMLEY 

R&  D  ASSOCIATES 

ATTN:  C  K  B  LEE 
ATTN:  D  SIMONS 

2  CYS  ATTN:  DOCUMENT  CONTROL 
ATTN:  FA  FIELD 

R  &  D  ASSOCIATES 

ATTN:  ADEVERILL 
ATTN:  C  KNOWLES 
ATTN:  J  THOMPSON 

RAND CORP 

ATTN:  P  DAVIS 
ATTN:  V  JACKSON 

RAND  CORP 

ATTN:  B BENNETT 

ROCKWELL  INTERNATIONAL  CORP 
ATTN:  J  HOWE 

S-CUBED 

ATTN  BPYATT 

SCIENCE  APPLICATIONS  INTL  CORP 
ATTN:  DOCUMENT  CONTROL 
ATTN:  ESWICK 
ATTN:  J  MARTIN 
ATTN  M  DRAKE 
ATTN  R  J  BEYSTER 

SCIENCE  APPLICATIONS  INTL  CORP 
ATTN  J  SHANNON 
ATTN  LGOURE 
ATTN  W LAYSON 
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SCIENCE  APPLICATIONS  INTL  CORP 
ATTN:  D  KAUL 

SCIENCE  APPLICATIONS  INTL  CORP 
ATTN:  RCRAVER 

TRW  INC 

ATTN:  DSCALLV 
ATTN:  R BURNETT 


